# We use gsynth 1.2.1, released on August 6, 2021.

# package
require("readr")
require("dplyr")
require("tidyr")
require("scales")
require("ggplot2")
require("gsynth")
require("panelView")
theme_set(theme_bw(base_size = 12))

# set working directory
setwd()

# load
utas.mrp <- read_csv(
  "MRP_PREFECTURE_UTAS.csv", locale = locale(encoding = "SJIS")
)

# recode
utas.mrp$ORDER <- rescale(utas.mrp$ORDER, to = c(1, 5), from = c(0, 1))

# time
utas.mrp <- mutate(utas.mrp, 
                   TIME = recode(YEAR, 
                                 `2003` = 1, `2007` = 2, `2009` = 3, `2012` = 4, 
                                 `2014` = 5, `2017` = 6, `2022` = 7, `2023` = 8))
table(utas.mrp$TIME)

# treatment variable
utas.mrp$TREATMENT <- 0
utas.mrp$TREATMENT[utas.mrp$YEAR == 2023 & utas.mrp$PREFECTURE %in% c(29, 35)] <- 1

# visualizing data
panelview(LDP ~ TREATMENT, data = utas.mrp, 
          index = c("PREFECTURE", "YEAR"), pre.post = TRUE)
panelview(LDP ~ TREATMENT, data = utas.mrp, 
          index = c("PREFECTURE", "YEAR"), type = "outcome")
panelview(ORDER ~ TREATMENT, data = utas.mrp, 
          index = c("PREFECTURE", "YEAR"), pre.post = TRUE)
panelview(ORDER ~ TREATMENT, data = utas.mrp, 
          index = c("PREFECTURE", "YEAR"), type = "outcome")

# generalized synthetic control
## feeling for LDP
## Nara
gsynth.nara1 <- gsynth(
  LDP ~ TREATMENT + DENSITY + COLLEGE_PERCENT + 
    FOREIGNER_PERCENT + PRIMERY_INDUSTRY_PERCENT + 
    SECONDARY_INDUSTRY_PERCENT + CRIME_RATE + 
    TURNOUT_PERCENT + LDP_PERCENT, 
  data = utas.mrp[utas.mrp$PREFECTURE != 35,],  
  index = c("PREFECTURE", "YEAR"), 
  se = TRUE, inference = "parametric", 
  parallel = TRUE, min.T0 = 7, 
  r = c(0, 5), CV = TRUE, force = "two-way", 
  nboots = 1000, seed = 12345
)
print(gsynth.nara1)
plot(gsynth.nara1, type = "gap", id = 29, main = "Nara")
gsynth.nara1.plot <- plot(
  gsynth.nara1, type = "counterfactual", id = 29, 
  main = "Nara", xlab = "Year", raw = "all") + 
  scale_color_manual(breaks = c("tr", "raw.co", "ct"), 
                     values = c("tr" = "black", 
                                "ct" = "black", 
                                "raw.co" = "gray90"), 
                     labels = c("tr" = "Treated", 
                                "ct" = "Estimated Y(0)", 
                                "raw.co" = "Controls")) + 
  scale_linetype_manual(breaks = c("tr", "raw.co", "ct"), 
                        values = c("tr" = "solid", 
                                   "ct" = "dashed", 
                                   "raw.co" = "solid"), 
                        labels = c("tr" = "Treated", 
                                   "ct" = "Estimated Y(0)", 
                                   "raw.co" = "Controls")) + 
  scale_size_manual(breaks = c("tr", "raw.co", "ct"), 
                    values = c("tr" = 0.6, "ct" = 0.6, "raw.co" = 0.5), 
                    labels = c("tr" = "Treated", 
                               "ct" = "Estimated Y(0)", 
                               "raw.co" = "Controls"))
plot(gsynth.nara1.plot)
## Yamaguchi
gsynth.yamaguchi1 <- gsynth(
  LDP ~ TREATMENT + DENSITY + COLLEGE_PERCENT + 
    FOREIGNER_PERCENT + PRIMERY_INDUSTRY_PERCENT + 
    SECONDARY_INDUSTRY_PERCENT + CRIME_RATE + 
    TURNOUT_PERCENT + LDP_PERCENT, 
  data = utas.mrp[utas.mrp$PREFECTURE != 29,],  
  index = c("PREFECTURE", "YEAR"), 
  se = TRUE, inference = "parametric", 
  parallel = TRUE, min.T0 = 7, 
  r = c(0, 5), CV = TRUE, force = "two-way", 
  nboots = 1000, seed = 12345
)
print(gsynth.yamaguchi1)
plot(gsynth.yamaguchi1, type = "gap", id = 35, main = "Yamaguchi")
gsynth.yamaguchi1.plot <- plot(
  gsynth.yamaguchi1, type = "counterfactual", id = 35, 
  main = "Yamaguchi", xlab = "Year", raw = "all") + 
  scale_color_manual(breaks = c("tr", "raw.co", "ct"), 
                     values = c("tr" = "black", 
                                "ct" = "black", 
                                "raw.co" = "gray90"), 
                     labels = c("tr" = "Treated", 
                                "ct" = "Estimated Y(0)", 
                                "raw.co" = "Controls")) + 
  scale_linetype_manual(breaks = c("tr", "raw.co", "ct"), 
                        values = c("tr" = "solid", 
                                   "ct" = "dashed", 
                                   "raw.co" = "solid"), 
                        labels = c("tr" = "Treated", 
                                   "ct" = "Estimated Y(0)", 
                                   "raw.co" = "Controls")) + 
  scale_size_manual(breaks = c("tr", "raw.co", "ct"), 
                    values = c("tr" = 0.6, "ct" = 0.6, "raw.co" = 0.5), 
                    labels = c("tr" = "Treated", 
                               "ct" = "Estimated Y(0)", 
                               "raw.co" = "Controls"))
plot(gsynth.yamaguchi1.plot)

## restricting rights to public safety
## Nara
gsynth.nara2 <- gsynth(
  ORDER ~ TREATMENT + DENSITY + COLLEGE_PERCENT + 
    FOREIGNER_PERCENT + PRIMERY_INDUSTRY_PERCENT + 
    SECONDARY_INDUSTRY_PERCENT + CRIME_RATE + 
    TURNOUT_PERCENT + LDP_PERCENT, 
  data = utas.mrp[utas.mrp$PREFECTURE != 35 & !is.na(utas.mrp$ORDER),],  
  index = c("PREFECTURE", "YEAR"), 
  se = TRUE, inference = "parametric", 
  parallel = TRUE, min.T0 = 6, 
  r = c(0, 4), CV = TRUE, force = "two-way", 
  nboots = 1000, seed = 12345
)
print(gsynth.nara2)
plot(gsynth.nara2, type = "gap", id = 29, main = "Nara")
gsynth.nara2.plot <- plot(
  gsynth.nara2, type = "counterfactual", id = 29, 
  main = "Nara", xlab = "Year", raw = "all") + 
  scale_color_manual(breaks = c("tr", "raw.co", "ct"), 
                     values = c("tr" = "black", 
                                "ct" = "black", 
                                "raw.co" = "gray90"), 
                     labels = c("tr" = "Treated", 
                                "ct" = "Estimated Y(0)", 
                                "raw.co" = "Controls")) + 
  scale_linetype_manual(breaks = c("tr", "raw.co", "ct"), 
                        values = c("tr" = "solid", 
                                   "ct" = "dashed", 
                                   "raw.co" = "solid"), 
                        labels = c("tr" = "Treated", 
                                   "ct" = "Estimated Y(0)", 
                                   "raw.co" = "Controls")) + 
  scale_size_manual(breaks = c("tr", "raw.co", "ct"), 
                    values = c("tr" = 0.6, "ct" = 0.6, "raw.co" = 0.5), 
                    labels = c("tr" = "Treated", 
                               "ct" = "Estimated Y(0)", 
                               "raw.co" = "Controls"))
plot(gsynth.nara2.plot)
## Yamaguchi
gsynth.yamaguchi2 <- gsynth(
  ORDER ~ TREATMENT + DENSITY + COLLEGE_PERCENT + 
    FOREIGNER_PERCENT + PRIMERY_INDUSTRY_PERCENT + 
    SECONDARY_INDUSTRY_PERCENT + CRIME_RATE + 
    TURNOUT_PERCENT + LDP_PERCENT, 
  data = utas.mrp[utas.mrp$PREFECTURE != 29 & !is.na(utas.mrp$ORDER),],  
  index = c("PREFECTURE", "YEAR"), 
  se = TRUE, inference = "parametric", 
  parallel = TRUE, min.T0 = 6, 
  r = c(0, 4), CV = TRUE, force = "two-way", 
  nboots = 1000, seed = 12345
)
print(gsynth.yamaguchi2)
plot(gsynth.yamaguchi2, type = "gap", id = 35, main = "Yamaguchi")
gsynth.yamaguchi2.plot <- plot(
  gsynth.yamaguchi2, type = "counterfactual", id = 35, 
  main = "Yamaguchi", xlab = "Year", raw = "all") + 
  scale_color_manual(breaks = c("tr", "raw.co", "ct"), 
                     values = c("tr" = "black", 
                                "ct" = "black", 
                                "raw.co" = "gray90"), 
                     labels = c("tr" = "Treated", 
                                "ct" = "Estimated Y(0)", 
                                "raw.co" = "Controls")) + 
  scale_linetype_manual(breaks = c("tr", "raw.co", "ct"), 
                        values = c("tr" = "solid", 
                                   "ct" = "dashed", 
                                   "raw.co" = "solid"), 
                        labels = c("tr" = "Treated", 
                                   "ct" = "Estimated Y(0)", 
                                   "raw.co" = "Controls")) + 
  scale_size_manual(breaks = c("tr", "raw.co", "ct"), 
                    values = c("tr" = 0.6, "ct" = 0.6, "raw.co" = 0.5), 
                    labels = c("tr" = "Treated", 
                               "ct" = "Estimated Y(0)", 
                               "raw.co" = "Controls"))
plot(gsynth.yamaguchi2.plot)
